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Three-component modeling of C-rich AGB star winds 

IV. Revised interpretation with improved numerical descriptions 
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ABSTRACT 

Models describing dust-driven winds are important for understanding the physical mecha- 
nism and properties of mass loss on the asymptotic giant branch. These models are becoming 
increasingly realistic with more detailed physics included, but also more computationally de- 
manding. The purpose of this study is to clarify to what extent the applied numerical approach 
affects resulting physical structures of modelled winds, and to discuss resulting changes. Fol- 
lowing the previously developed radiation hydrodynamic model - which includes descriptions 
for time-dependent dust formation and gas-dust drift - and using its physical assumptions and 
parameters, numerical improvements are introduced. Impacts of the so-called adaptive grid 
equation and advection schemes are assessed from models calculated with different numeri- 
cal setups. Results show that wind models are strongly influenced by numerical imprecision, 
displaying differences in calculated physical properties of up to one hundred per cent. Using a 
non-adaptive grid, models become periodic (in multiples of stellar pulsation periods), instead 
of irregular, as obtained previously. Furthermore, the numerical improvements reveal changes 
in physical structures. The influence of gas-dust drift is confirmed to be highly important, in 
particular for the dust component. Gas and dust are less tightly coupled than previously, and 
drastically larger amounts of dust form assuming drift. 
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1 INTRODUCTION 

Dust-driven winds off of asymptotic giant branch (AGB) stars are 
complex physical phenomena that form in a highly non-linear in- 
teraction between a strong radiation field, a cool gas and condi- 
tions suitable for efficient dust formation. It is the radiation pres- 
sure on dust grains that with the support of levitated regions of 
an enclosed pulsating st ar is believed t o form and drive more mas- 
sive winds. Starting with lBowenI 1 119881) a variety of time-dependent 
wind models, of increasing complexity, have been created in order 
to understand AGB star winds physically. An overview of proper- 
ties an d capabilities of recent w ind mode l s is gi ven by, e.g.. lHofnen 
( I2OO5L and references therein). IWillsonI ( 120000 reviews mass loss 
from cool stars in more general terms. 

An important feature of self-consistent wind models, as pre- 
sented here, is an adaptive grid equation. This equation distributes 
a fixed number of gridpoints across the model domain according 
to thos e physical properties that are required to be spatially re- 
solved l lDorfi&Drun^l 19871 henceforth DD87). The intent to use 
an adaptive grid is sound because shocks form in radially steep den- 
sity gradients of long period variables. Shocks that assist in pro- 
viding appropriate conditions for dust formation and hence wind 
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acceleration. It has earlier been found that a vast majority of all 
wind models using such a grid e quation have formed i rregular 
outflow patterns; for exam ples see: iHofner & Dorfil l ll997l fig. 1). 
ISandin & HofneJ (12003 j, henceforth P aper I, fig. 2). iHofner et"ai] 
( 120031 fig. 2).ISandiii & HofneJ ( l2003bl henceforth Paper II, fig. 5), 
and iMattsson et al.l 1 120071 figs. 5-7). Such irregularities could be 
understood as a result of physically (mildly) chaotic interactions, 
but they might also be a consequence of insufficient numerical ac- 
curacy. 

Few stu dies addres s the u se of the adaptive grid equation. To 
mention two iKiirschneii ( 119941) compares results between models 
using first and second order accurate advection schemes; conclud- 
ing that a second order advection is important. iDorfi et al.l ( l200q . 
henceforth D06) introduce both an improved discretization scheme 
and an improved advection scheme. In addition to the advection and 
discretization schemes, other "numerical properties" also have the 
potential of affecting the physical outcome of models. Such exam- 
ples are how the adaptive grid equation is used, and the adopted grid 
resolution. There is no thorough study of wind models available in 
the literature, where the influence on results due to the numerical 
approach is explored. For models using the adaptive grid equation 
there is no such study available at all. This is an important issue of 
concern, since a model interpretation is uncertain when it is unclear 
to what degree the outcome is affected by numerical inaccuracy. 
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This paper is based on the work presente d in the previous 
papers in the series, i.e., Paper I, Paper II, and I Sandin & Hofnen 
(2004. henceforth Paper III); a comprehensive summary is given in 
[S^andin (2003). The purpose of this paper is twofold. Firstly to fill 
the gap of a numerical study by exploring and determining a mod- 
elling approach by which numerical effects on results are separated 
and minimized. All physical parameters are thus fixed to previously 
used values. Secondly, the improved numerical approach leads to 
significant changes in the outcome, which are new and well worth 
addressing. For instance, when the adaptive grid does not track 
shocks outflow structures turn out to be periodic, instead of irreg- 
ular as before. Assuming drift the wind also becomes more tenu- 
ous, simultaneously drastically larger amounts of dust are formed. 
Issues of the numerical method and physical results are treated to- 
gether in order to emphasise how strongly they depend on each 
other. Results should, finally, also be valuable when improving re- 
lated studies working with mas s loss rates and yields of dust from 
AGB stars further, such as, e.g.. lFerrarotti & Gaill ( 120060 . 

Numerical modifications and improvements are introduced, 
after first summarising the physical and numerical setup of the cur- 
rent AGB star wind model, in Sect.[2] The modelling procedure and 
model setup are then described in Sect. (3) where results and code 
consistency checks are also presented. The improved numerics in 
several cases gives rise to drastic changes to the physical structure. 
These changes and the revised outcome are discussed, after differ- 
ences of the numerical modifications are treated, in Sect. |4] The 
paper is closed with conclusions in Sect.|5] 



2 MODEL FEATURES AND NUMERICAL 
IMPROVEMENTS 

Following the three previous papers in this series, a distinction is 
made between three interacting physical components in the wind; 
gas, dust, and a radiation field. The system is described by coupled 
conservation equations, which account for exchange of mass, en- 
ergy, and momentum between all three components - forming the 
RHD3-system of equations (Radiation Hydrodynamics Dust and 
Drift). A thorough description of the physical system, the gas-dust 
interaction, and the numerical method was given in Paper I (and 
references therein; where also a constant gas opacity was used), 
while Paper II added effects of stellar pulsations and grey gas opac- 
ities. Paper III, finally, introduced effects of gas-dust drift in the 
dust formation process. The work presented here is based on this, 
most recent, physical formulation. All physical assumptions and 
parameters, as well as the general aspects of the numerical method, 
are unchanged with respect to these three papers, in order to fo- 
cus on effects due to the numerical modifications studied here. In 
comparison to other AGB star wind models, these are based on - 
and have most facto rs in common with - the models presented by 
iHofneret all d 19981) . 

Gas-dust drift is considered ubiquitous in this article. How- 
ever, since no other known AGB star wind model can handle drift 
without serious limitations, more restrictive non-drift (position cou- 
pled [PC]) models are also treated. This approach allows a more 
precise study than before of differences, both due to the numerical 
approach, and between drift and PC models. 

New to this paper are modified descriptions of the spatial dis- 
cretization and advection schemes, and how the adaptive grid equa- 
tion is used. These modifications are introduced in Sects. [2^2.41 
after a summary of the numerical method used so far is presented in 
Sect. 12. II Important issues related to the size of the model domain 



Table 1. Glossary of used abbreviations and symbols 



Abbreviations: 



term 



description 



AGB asymptotic giant branch 

RHD radiation hydrodynamics 

RHDD radiation hydrodynamics & time dependent dust formation 

RHD3 radiation hydrodynamics & time dependent dust formation 

- including drift 
PC position coupling (non-drift) 

Advection schemes; cf. Sect. |2.3| 
vL van Leer (second order) 
VWvL volume weighted van Leer (second order) 
PPM piecewise parabolic method (third order) 

Grid types - use of the adaptive grid equation 
A adaptive and logarithmic distribution 
L non-adaptive and logarithmic distribution 
U non-adaptive and uniform distribution 

Symbols and properties: 



symbol 


unit 


description 


r 


cm 


radius 


P 


gcm-3 


gas density 


e 


ergg-l 


specitic internal energy of the gas 


Tg 


K 


gas temperature 


Kg 


cm2g-l 


gas opacity 


u 


cms~^ 


gas velocity 


J* 


s~^cm^^ 


net grain formation rate per volume 


rn' 


s-i 


net grain growth rate 


Pi 


gcm-3 


dust density 


VD 


cms~l 


drift velocity (size averaged) 


Ki 


cm-3 


moments of the grain size 
distribution; ^ j ^ 3 


/cond 




degree of condensation 


Pi/p 




dust-to-gas density ratio 


I'd 


cm 


mean grain radius 


R 




adaptive grid equation resolution function 


w 




adaptive grid equation weights for R 


Model 


input parameters 


L* 


ergs-l 


stellar luminosity 


M, 


g 


stellar mass 


fl* 


cm 


stellar radius 


Teff 


K 


effective temperature 


=c/eo 




elemental abundance of C relative to O 


P 


s 


stellar pulsations: piston period 


Awp 


cms"-"^ 


stellar pulsations: piston amplitude 


^int 


cm 


radial location of the inner boundary 


fcxt 


cm 


radial location of the outer boundary 


rig 




number of gridpoints 


Properties only calculated at the outer boundaiy 


M 


M0yr-i 


mass loss rate 


Uoo 


cms~l 


terminal velocity (final wind velocity) 


f 




relative fluctuation amplitude 


{q) 




temporal mean of quantity q 


CTs 




standard deviation 



and grid resolution are additionally addressed in Sect. 12. 51 All sym- 
bols and abbreviations used in this paper at more than one location 
are for a quick reference summarised in Table[T] 
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2.1 Numerical method - up to now 

For a detailed description of the numerical method used so far see 
Paper I, sect. 3.1, and references therein. Those features that are 
modified in this article are summarised in the following. 

An adaptive grid equation iDorfi & Drurv 1987) distributes 
gridpoints on the model domain by resolving gradients of selected 
quantities. Hence, gridpoints move back and forth according to the 
definition of a so-called grid resolution function R (cf. ibid., eq. 3). 
The discretized and dimensionless form of R can be written. 



R^^ 



\ 



M 



X, h 



h 



F,. 



ri+i 



(1) 



where Xi is a spatial scale length at gridpoint i, r is the radius, 
Fj,i the scaling factor for the physical quantity fj,i to be resolved 
(totally there are A'l such quantities for each gridpoint), and Wj are 
grid weights normally set to 1 .0 for the resolved quantities in stellar 
wind models (these grid weights are not written out explicitly in 
the original expression). Like before R is set to be determined by 
the thermal (internal) energy (e) and the gas density (p). Abrupt 
changes in the distribution of gridpoints are, moreover, controlled 
using two smoothing factors. The temporal smoothing factor is set 
to Tg — 10^ s, which is smaller than, e.g., dynamical or dust related 
time scales in the wind, allowing the grid to freely adapt to physical 
features. The spatial smoothing factor is set to q = 2. Typically the 
gas density drops by a factor of 1.0-1.5 dex across a stronger shock 
why at least 10-15 gridpoints are required to achieve a sufficient 
local refinement (with a = 2, cf. DD87, p. 182). 

Effects of stellar pulsations on the atmosphere and wind are 
described using a sinusoidal, radially varying inner boundary (of 
period P; i.e., a piston), located at about ri„t = 0.91 -R*; above the 
region where the K-mechanism supposedly originates . An inflow 
of mas s through the i nner boundary - like used by, e.g.. lSimis et al.l 
( I2OOIL sect. 4.2) and lWoitkd ( l200d henceforth W06, sect. 2.5.1) 
- is as before not permitted. The typical amount of mass in the 
modelled envelope is about 0.002Mq. Since the model domain 
is rather quickly depleted of material due to the wind, long-term 
modelling, covering thousands of years, is (currently) not possible. 

The system of equations is discretized in a volume-integrated 
conservation form on a staggered mesh. Advection of mass, energy, 
and momentum between grid cells is de scribed using a s econd or- 
der accurate discretization according to Ivan Leen 1 119771) . To keep 
the formulation fully implicit the temporal discretization, of all 
terms in all equations, is always first order accurate. The equations 
are solved implicitly using a Newton-Raphson algorithm where the 
Jacobian of the system is inverted by the Henyey method. 

2.2 Discretization schemes - introducing weighted means 

The equations of radiation hydrodynamics in spherical geometry 
are discretized on a staggered mesh, i.e., scalar quantities such as 
the density are centred within a computational cell. Vector quanti- 
ties, such as the velocity and radius, are localised to cell boundaries. 
When combining scalar and vector type quantities one must adopt 
a discretization scheme, which has up to now always been an arith- 
metic average of variables at neighbouring gridpoints (regardless 
of geometry). E.g., 



= 2 (''^ + '■^+1) 



and p, = -(p,+ i+p. 



(2) 



where r^ is the radius on the cell boundary at gridpoint i, and r^_^i 
the averaged radius at the cell centre (ri+i < r^). Other quantities 



are treated analogously. As pointed out by, e.g., D06 (see sect. 2.1) 
a more consistent discretization scheme for a spherical geometry is 
achieved if quantities instead are weighted over the volume V of 
each grid cell. Thus, the relations in Eq. (|2j are replaced with. 



r,+ 1 = V 2 I'-J 



+ rf+i) and 



pi = 



^(p,+ iK+i+P,_iK_i) 



(3) 



This improved discretization scheme has been implemented in 
all terms, and is used with all models using the volume weighted 
van Leer advection scheme (see next subsection). However, unlike 
the choice of the advection scheme (or use of the adaptive grid; 
cf. Sect. 12.4b . the choice of the discretization scheme has in tests 
carried out for this paper been found to have a negligible influence 
on the physical structure of wind models. 



2.3 Advection schemes - increasing the accuracy 

The adopted advection, describing transport of matter, energy, and 
momentum between grid cells, has up to now b een the spatially 
second order accurate (van Leer [vL]) scheme of Ivan Leeii i ll 9771 . 
eq. 671J. In order to study the influence on stellar winds of higher 
accuracy in the advection, both the volume weighted second or- 
der (VWvL) scheme of D06 as well as th e third order piecewise 
parabolic method of IColella & Woodward 1 1984 , PPM, which is 
also volume weighted; in the version of Iwinkler & Normanll 19861 
- henceforth WN86 - sect. V-Co have been implemented. A tem- 
poral discretization only first order accurate motivates this use of 
higher order spatial discretization in the advection scheme, to reach 
an overall better accuracy. 

Using the (integrated) mass in the advection term discretiza- 
tion (see, e.g., WN86, sect. V-B), the system Jacobian (cf., e.g., 
[Dorfi 1998, sect. 5.3) becomes block penta-, hexa-, and nona- 
diagonal for the vL, VWvL, and PPM schemes, respectively. For 
the PPM scheme variables at the neighbouring gridpoints (i — 4, i — 
3, . . . , i, . . . , J + 4) enter the discrete equations. Compared to the 
vL scheme the time penalties for each Gaussian elimination of the 
Jacobian are 20 per cent (VWvL) and 80 per cent (PPM). 

In difference to the gas component it does not make sense to 
evaluate an integrated mass for the dust, which is why drift models 
always must be advected in the dust component using a velocity. 
With a velocity instead of a mass in the advection term discretiza- 
tion the corresponding Jacobian is block hexa-, hepta- and nona- 
diagonal. Time penalties are 20 per cent (vL), 40 per cent (VWvL), 
and 80 per cent (PPM), respectively; when compared to the mass- 
advected vL scheme. It should be noted, however, that since the 
time-dependent iteration history changes with different advection 
schemes it is difficult to say exactly how the total model compu- 
tation time changes. That the VWvL and PPM schemes have been 
correctly implemented is seen in Sect. 14.1.71 where both give very 
similar structures, although PPM achieves slightly more detail. 

Up to now it has always been necessary to run adaptive grid 
wind models using the mass advection formulation in the gas equa- 
tions, or the same equations fail to converge properly (resulting in 



^ The drift models in Paper I use a spatially first order accurate advection 
scheme, referred to as Donor Cell. 

^ Which inclusion is a capital task due to the large number of analytical 
derivatives required in an implicit formulation. 
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very small time steps). In this context D06 (in a different appli- 
cation) successfully switch from advection by mass to advection 
by velocity, after introducing the volume weighted discretization 
(Sect.|2J2j and VWvL advection scheme. The same procedure does 
not work (in general) with wind models using an adaptive grid. 
However, with a non-adaptive grid (see next subsection) wind mod- 
els can too switch to the velocity advection formulation. Compar- 
isons between mass and velocity advected models show negligible 
differences, which is why the mass formulation is used to save time 
in the calculation of PC models. 



2.4 The use of an adaptive or non-adaptive grid equation 

Narrow shocks form in the steep density gradients of a pulsating 
AGB star atmosphere. The use of adaptive mesh refinement or an 
adaptive grid equation is motivated in order to resolve the phys- 
ical structure of such shocks. There are, however, several draw- 
backs with possibly major consequences for the physical structure 
of the modelled region when an adaptive grid equation is used. 
Such drawbacks have not been addressed in any detail before and 
are treated in the following subsection. An approach using a non- 
adaptive grid (equation) is thereafter discussed in Sect. 12.4.2] 



(where i?, is the stellar radius). Incidentally such numerical arte- 
facts always occur in front of shocks where there is very little dust 
(and few gridpoints). If and when the dust drifts ahead of the gas 
these numerical artefacts together with a locally decreased resolu- 
tion can cause major problems to the overall physical structure of 
the wind, cf. Sect. 14.2.11 This numerical issue makes drift mod- 
els difficult to model with an adaptive grid equation. With a non- 
adaptive grid, which does not track shocks, the advection velocity 
is non-zero in regions with dust - with the possible exception of 
cool inner regions - and this problem is greatly reduced (artificial 
dust mass diffusion is still required occasionally). 

A wind model covers a large spatial domain, so far typically 
0.91-25 i?*, with plenty of structures in both the gas and dust, 
which can all be important to the overall structure of the accel- 
erating wind. In particular drift models form more structures in the 
dust than PC models do (cf., e.g., sect. 5.2 in Paper II). A large 
number of structures, which cannot all be resolved simultaneously, 
and the difficult numerical conditions for modelling drift are both 
good motives for using a non-adaptive grid. Still, one strong argu- 
ment to keep the grid equation (in a non-adaptive setup) is that it 
allows a simple radial displacement of the inner boundary to simu- 
late pulsations (cf., e.g.. Paper II, sect. 2.3). 



2.4.1 Issues of concern with an adaptive grid equation 

While certain regions are adequately resolved through an appro- 
priate definition of the grid resolution function, unresolved regions 
may, on the contrary, be "gridpoint exhausted". In particular if the 
model domain covers a spatially large region with a large number 
of structures. That is, structures may fail to form in a region where 
they could have formed if the same region was better resolved. This 
is of particular concern in drift models where large variations form 
in dust quantities between regions resolved by the grid equation. 
There are two immediate solutions to this problem, adding more 
gridpoints or using a different resolution function, but these solu- 
tions do not work well. Choosing a different resolution function in 
order to resolve quantities of the dust is difficult since dust prop- 
erties primarily vary on shorter time scales than those of the gas; 
as the dust component lacks a pressure component and reacts very 
quickly to changes in the radiation field. Resolving both gas and 
dust features would at times require most of the model domain to 
be resolved simultaneously, disapproving its use. 

Adding significantly more gridpoints is in general not an op- 
tion, as the relative difference in the amount of mass between grid 
cells may become too large (since most gridpoints will be resolving 
shocks further), introducing new numerical problems. However, 
models using a mildly increased(/decreased) number of gridpoints 
should reproduce average properties, if the numerical scheme is 
stable (cf. Sect. |2.5l l. A third option is to use smaller grid weights 
in the grid resolution function (wj in Eq.lTJ, cf. Sect. l4.1.2l 

There is a critical issue with the grid equation that is en- 
tirely limited to drift models. In a grid cell with a small amount 
of mass (like in the dust component) there is the possibility that the 
amount of numerical diffusion due to advection becomes too small 
to prevent numerical errors from growing out of control (see, e.g., 
sect. IV-C in WN86, and sect. 3.2 in Paper I). This is the case in 
regions where gridpoints are moving in both directions to resolve 
physical (or numerical) features elsewhere and the local relative 
advection velocity becomes close to zero. In order to remedy the 
problem some artificial mass diffusion can be added to the dust 
component. An illustration of this issue can be seen in the dust ve- 
locities shown in Paper I, figs. 2b & j at r = 15 -R* and r — 25 Ri, 



2.4.2 Introducing a non-adaptive grid equation 

On top of a slightly moving grid due to a displaced inner boundary 
(see above) the basic distribution of gridpoints can (most easily) be 
made logarithmic or uniform (henceforth grid types L and U). A 
logarithmic grid has the advantage that narrow shocks in the stellar 
centre are better resolved than with a uniform grid, also with fewer 
gridpoints. A logarithmic grid is therefore preferred over a uniform; 
the latter is, however, used for comparison. 

Two relevant properties of the adaptive grid equation in this 
context are the resolution function R (Eq.[TJ and the gridpoint con- 
centration n (see DD87, eq. 2). The grid equation is non-adaptive 
- and not tracking shocks - if _R is set to 1. n is defined by, 



Xi 



Ti ■ 



■ ri+i 



and Xl , 



X\] , 



c, 



(4) 



where ^ a: ^ 1. A logarithmic grid is achieved if the scale length 
is set to X = Xl (a; = 1). Similarly a uniform distribution results if 
the scale length instead is set to X = Xu, where C is a constant. If 
X = Xl and a; = 0.5 the grid will be very nearly uniform and non- 
adaptive (also with R ^ 1), since the nominator and denominator 
of n differ by orders of magnitude; using typical values of r. This 
setting is used with models using a spatially uniform grid. 

With a non-adaptive grid time steps will necessarily be shorter 
on average since the "stationary" temporal periods found in models 
using an adaptive grid vanish. It simultaneously turns out that the 
required number of iterations in each time step of the non-adaptive 
model on the average is lower. In total the computational time re- 
quired to reach a specific evolved "age" of a model is slightly higher 
for a non-adaptive model compared to the adaptive counterpart. As 
is found in Sect. 13.31 all presented models using a non-adaptive 
grid form outflows with a periodic character instead of irregular, 
as was found previously. It is thus not necessary anymore to evolve 
a model to reach a high age in order to get sensibly temporally av- 
eraged quantities. 
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2.5 Size of model domain and grid resolution 

Models typically span a region from about rint = 0.91 -R* to 
rcxi ~ 25-60 i?*, covering the atmosphere and wind acceleration 
region. In earlier papers in the series, the outer boundary was fixed 
at about r^xt = 25 i?*. W06 in his study places the outer boundary 
at Text = 10 i?*. The argument behind these values is that the wind 
has reached its terminal velocity within these domains. 

The default number of gridpoints is Ug — 500 for both models 
using an adaptive and non-adaptive grid; this is also the same num- 
ber used previously. An issue related to the grid resolution is the 
amount of added artificial tensor viscosity. This amount is defined 
using a length scale (cf., e.g., Paper I, eqs. 13 & 14), 



^relative = rf, VS. /flxcd — Tq f , 



(5) 



where / is a constant that defines the shock width and ro another 
constant. In order to always resolve shocks in the innermost re- 
gion with at least one gridpoint (in models using a logarithmic non- 
adaptive grid, with rig — 500 & Text ~ 50 /?*) /relative is used with 
/ = 7.0 X 10~^. In models using a uniform grid ioxcd is used in- 
stead, with ro = 5 7?* . Differences in the physical structure due to 
the chosen value on the length scale are found to be negligible for 
/ < 2.0 X 10^^ (using P13C16U6„„=7oo). With larger values the 
physical structure changes from periodic to irregular 

Finally, the influence of the adopted number of gridpoints 
on the physical structure has not been studied in any detail with 
these stellar wind models earlier. Ideally models using an increased 
number of gridpoints should become more alike, if the numeri- 
cal method (and physical circumstances) permits. In addition it 
is useful to see how dependent the physical structure and tempo- 
rally averaged quantities are on a decreasing number of gridpoints. 
An accurate spatial discretization (Sect.|2J2]l and advection scheme 
(Sect. |23l l are of increased value with less dense grids. To find out 
one set of model parameters has been run using models with four 
different number of gridpoints, cf. Sect. l3.3l 



3 MODELLING PROCEDURE AND RESULTS 

In this section the modelling procedure is, at first, reviewed in 
Sect. 13.11 followed by a description of the chosen physical setup 
and selected model parameters in Sect. 13.21 Results are then pre- 
sented in Sect. 13.31 A discussion of models using a constant gas 
opacity is thereafter given in Sect. 13.41 The outcome of some tests 
for code validation is provided in Sect. 13. 51 



3.1 Modelling procedure 

The modelling procedure is as follows. A wind model is started 
from a hydrostatic dust-free initial model, where the outer bound- 
ary is located at about rcxi = 2 i?*. All dust equations and terms are 
then switched on simultaneously. Dust starts to form in the outer 
cooler domain, initiating an outward motion of dust and gas. The 
expansion is followed by the grid to a pre-defined radius (of about 
Text = 25-60 -R*), where the outer boundary is fixed and subsequent 
outflow allowed. Thereafter the actual wind is modelled. 

A wind model is evolved for a time interval of about 
50-200 P. The modelled time interval cannot be much shorter than 
50 P (when Text = 25 /?*), which is the typical time needed for tran- 
sients in the expansion phase to leave the model domain through 
the outer boundary. Temporally averaged quantities should be mea- 
sured only after this time; and larger radial domains consequently 



require longer modelled time intervals. Measuring averaged quan- 
tities extended time intervals are unimportant for models using a 
non-adaptive grid, since these in all cases studied here are found to 
develop periodic outflows. For such models precise averaged quan- 
tities can be calculated over an interval of a few piston periods. The 
situation is different with models using an adaptive grid, which for 
most computed cases develop an irregular structure, like before. 
Such structures require a longer time interval to sensibly measure 
averaged quantities. 

Models using a Planck mean gas opacity (see below) all de- 
velop a 'low mass' envelope containing about 10~^ ^Q- ^^^ '■™^ 
interval used to compute averaged quantities must be short enough 
that the model domain is not depleted of mass (also see sect. 4.2 
in Paper II). The intervals studied in Sect. l3.3l are selected at times 
before a significant fraction of the envelope is lost; in comparison 
to a model with periodic structures it is in general more difficult to 
define a suitable short time interval with irregular structures. 



3.2 Physics setup and selection of model parameters 

All models in this paper are selected with the purpose of studying 
relative changes compared to previously calculated wind models. 
This approach allows a quantitative and qualitative estimate of the 
importance of the adopted numerical modifications. 

Highly critical to properties such as the developed density 
structure is the gas opacity used. Models adopting Planck mean 
absorption coefficients of molecular data (Planck mean opaci- 
ties/models) on the one hand result in much more realistic density 
structures, compared to models calculated with a constant gas opac- 
ity (constant-opacity models; see, e.g., sect. 2.2 in Paper II). On 
the other hand Planck mean models are less realistic than models 
which use frequenc y-dependent opacities (for both the gas and dust, 
iHofner et alj2003r) . Although frequency-dependent opacities intro- 
duce another level of complexity to the problem, and are also very 
computationally demanding. In the current context Planck mean 
gas opacities are used to avoid these problems. 

In comparison to Paper III, where the influence of various lev- 
els of drift inclusion in the dust formation equations was studied, all 
drift models here use the complete drift-dependent dust formation 
description (marked with the suffix "—dvs" in ibid.). One limita- 
tion of this formulation is that dust formation in principle is grain 
size dependent through the drift velocity (cf. Sect. |4.2.4l . The com- 
putational requirements of a study using a grain size distribution are 
much larger than with a mean grain radius, which is why such an 
approach is out of scope of this paper. Still, no other known time- 
dependent AGB star wind model includes drift-dependent dust for- 
mation, using even a mean grain size. All remaining physical pa- 
rameters of the gas and dust are unchanged. 

The physical equations solved are highly non-linear, as is the 
outcome. A larger number of models have accordingly been cal- 
culated to provide better statistics. In order to enable a straightfor- 
ward comparison with earlier results, the same model parameters 
used in Paper II (and Paper III) have also been used here; see Ta- 
ble |2] Models for which a wind did not form earlier (PC models 
with a terminal velocity Uao < 10 kms^^) are also left out here, 
since they are not expected to form a wind now either. One model, 
P13C16U6, is used here in a thorough study of the influence of the 
number of gridpoints. 
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Table 2. Model parameters, see Sect. |3.2| for further details. The name of 
the model is given in Col. 1. The following five columns specify: the stel- 
lar luminosity L^-, effective temperature T^ff, pulsation period P, pulsation 
amplitude A «p, and carbon-to-oxygen ratio ^c/sq. The last column shows 
the fraction of the stellar mass contained in the radial domain of the initial 
model. The stellar mass M-,, is set to 1.0 A/q in all models. These sets of 
model parameters have been taken from table 2 in Paper II. 



model 




Teff 

[K] 


P 

[d] 


A Up 
[km/s] 


^c/eo 


A/e 

[%] 


P10C16U6 


1.0 


X 10* 


2790 


525 


6 


1.6 


0.16 


P10C18U4 


1.0 


X lO* 


2790 


525 


4 


1.8 


0.13 


P10C18U6 


1.0 


X 10^ 


2790 


525 


6 


1.8 


0.13 


P13C14U6 


1.3 


X 10-1 


2700 


650 


6 


1.4 


0.19 


P13C16U4 


1.3 


X 10" 


2700 


650 


4 


1.6 


0.16 


P13C16U6 


1.3 


X 10* 


2700 


650 


6 


1.6 


0.16 



3.3 Results 

In order to assess the importance of the numerical modifications, 
each set of model parameters has been used in a number of differ- 
ent setups. Both drift and PC models have been calculated for all 
sets. All models have additionally been calculated using at least two 
out of three advection schemes - marked with either vL, VWvL, or 
PPM (see Sect. |2.3l l. Models using an adaptive grid (grid type A) 
have been recomputed for PC models using vL and VWvL advec- 
tion. Both to see what the influence of the advection scheme is, and 
to see how new values of recomputed vL-models match previously 
calculated values - as a consistency check (cf. Sect. 13. 5b . Drift mod- 
els have not been calculated at all here with an adaptive grid due to 
the numerical difficulties mentioned in Sect. I2.4.TI Moreover, all 
models using a non-adaptive grid have been calculated using both 
a logarithmic (L) and a uniform (U) grid type. Since the emphasis 
has been to resolve shocks, only results of models using a logarith- 
mic grid are presented; except for P13C16U6. 

The wind is characterised through temporally averaged prop- 
erties, such as the mass loss rate (A/) and terminal velocity (woo) 
for the gas. These two quantities together give the gas density (p). 
The dust is characterised through four properties: the degree of con- 
densation (/cond), dust-to-gas density ratio (pj/p), mean grain radius 
(rd), and drift velocity (vd) (only for drift models). The degree of 
condensation is defined as. 



f -PL 

/cond — ,„, 



K. 



K^ + nc' 



(6) 



where pi = miK-j, is the dust density, and K^ the third moment 
of the grain size distribution (mi is the dust grain monomer mass; 
cf., e.g.. Paper I, sect. 2.1); and p[P' the total density of condensible 
matter (present in both the gas and dust phases); n^ the total number 
density of condensible material in the gas phase. Providing a mea- 
sure for the variability of the model structure, each outflow prop- 
erty iq) is attended by a relative fluctuation amplitude (f = (Ts/g; 
where CTs is the [sample] standard deviation of the quantity q in the 
measured time interval). 

Results for models using the selected model parameters, dif- 
ferent advection schemes, and grid types are presented in Table[3] 
In order to emphasise the increased accuracy of models using a 
non-adaptive grid, the conesponding values are written out with 
three significant digits. For each set of model parameters PC mod- 
els are listed before drift models. Additionally, for each set the two 
drift and PC models calculated using the best available advection 



scheme and a non-adaptive grid are treated and marked as refer- 
ence models (REF.). Other models in a set are marked in boldface 
if a property deviates by more than 10 per cent (i.e., significantly) 
from the corresponding reference drift model value. Likewise, PC 
models are underlined if a property deviates by more than 10 per 
cent from the respective PC reference model value. 

In order to clarify differences all properties are in addition to 
Table[3]also illustrated in Fig.[5j for absolute values (top panels) as 
well as values normalised to the respective value of each PC model 
(bottom panels). For each quantity (g) error bars in the figure show 
the property. 



= 10^ X logio [fil) 



-riM)] ■ 



(7) 



where fdi (g) indicates a normalisation to the respective drift model 
value. Due to significant variations of the relative fluctuation ampli- 
tude for separate properties in different models this representation 
allows a visually better comparison. 

The influence of the grid resolution is studied for model 
P13C16U6 using four different numbers of gridpoints: ng=100, 
300, 500, and 700 (the outer boundary is fixed at Text ^ 50 i?*). 
This model was part of a detailed study in Paper III. It is chosen 
for discerning differences arising from different resolutions, since 
it shows fairly small differences due to different advection schemes 
and grid types with ng=500. Results are given in Table|4l The num- 
ber of gridpoints is printed in subscript following the model name 
when referring to these models. 



3.4 Comparing outcome with (grey) constant-opacity models 

Several constant-opacity models have also been calculated with all 
numerical modifications applied. As a contrast to the Planck mean 
models of this study a majority of the constant-opacity models de- 
velop an iiTegular structure, both with an adaptive and non-adaptive 
grid. Differences between outflow properties of the two models ap- 
pear to be small. That these models form irregular structures, also 
when using a non-adaptive grid, is in good agreement with the lat- 
est results of W06. He computes constant-opacity one-dimensional 
wind models using a code with adaptive mesh refinement. 

Moreover, by calculating the model P13C16U6 with a pre- 
defined constant gas opacity - such as Kg = 4 x 10^^ g~^cm^ - 
an envelope results that has a very similar amount of mass as the 
original model. Compared to the periodic structure of the original 
Planck mean model this constant-opacity model forms a wind with 
irregular structures. The conclusion is that constant-opacity mod- 
els are more "chaotic" and less sensitive to numerical details than 
Planck mean models are; because of the fewer degrees of freedom, 
possibly in combination with very massive envelopes. 



3.5 Consistency checks - code validation 

Large modifications have been introduced to the code since its last 
application in Paper III. It is, hence, important to compare new val- 
ues with old to see how similar they are. New models have been 
calculated for all sets of PC model parameters using an adaptive 
grid; these are marked with vL-A in Table [3] (Cols. 3 & 4). The 
values should be compared with those values of the correspond- 
ing model in Paper II (see table 5 therein). For the PlO-prefixed 
models the agreement is within 10 per cent for all three models and 
quantities, except for the mass loss rate of P10C18U6, which is 
25 per cent larger here. The disagreement is bigger for the P13- 
prefixed models. P13C16U6, first, shows the same agreement as 
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Table 3. Temporally averaged quantities at the outer boundary; see Sect. 13.31 Several models with different numerical setups are presented for each set of 
modelled parameters. From the left the first five columns specify input properties: model name, PC or drift (p/d), advection scheme (A), grid type (g), and 
radius at the outer boundary (rext). Remaining columns characterise the outflow in seven properties: the mass loss rate (M), terminal velocity (mcxd). degree 
of condensation (/cond)i dust/gas density ratio (Pci/p), dust radius (rj), drift velocity (tiu) (only for drift models), and type of outflow (type) for the gas:dust. 
Advection schemes (vL, VWvL and PPM) are described in Sect. 12.31 and grid types (A, L and U) likewise in Sect. 12.41 The types of outflow structures that 
form are classified as: s, stationary: i, irregular; I p, periodic; and I q quasi-periodic. I (S N) indicates the (multi-)periodicity of the gas/dust outflow in the unit 
of the piston period P. Models using an adaptive grid all show the same type of structure in both the gas and dust. In addition a relative fluctuation amplitude 
f (= <Js/q; (Ts is the standard deviation) is specified for each quantity q (cf. Sect. 13.31 . The most accurately calculated drift and PC (reference) models are 
marked with REF. in Col. 1. Values shown in boldface/underline indicate a significant difference (Jj 10 per cent) from the conesponding diift/PC reference 
model value. 



model 


p/d 


A 


9 


r-ext 


lOS 

[Mq 


yr-l] 
lO'V 


{v 
[km 


s-1] 
lO^r 


(/cond) 

[%] 

lO'V 


(Pd/p) 
[10-4] 
lO^r 


102{rd) 
[^lm\ 
lO^r 


[kms-i] 


type 


P10C16U6 


PC 


vL 


A 


25 


M 


480 


13 


85 


25 


160 


M 


160 


3.0 


150 






i 






VWvL 


A 


25 


M 


760 


M 


93 


28 


170 


M 


180 


M 


180 






i 






vL 


L 


32 


2.45 


7.1 


10.4 


1.2 


18.2 


6.6 


6.16 


6.1 


3.02 


5.0 






s:s 


REF. 




VWvL 


L 


32 


1.69 


21 


10.9 


2.5 


18.7 


2.7 


6.42 


2.7 


2.79 


1.8 






s:q 




drift 


vL 


L 


40 


1.31 


440 


14.5 


48 


27.7 


990 


11.4 


1200 


3.30 


460 


8.82 


460 


5q:5q 


REF. 




VWvL 


L 


40 


1.24 
1.2 


270 
470 


15.8 
15 


16 

41 


38.7 
17 


760 
140 


16.4 
7.7 


1000 
140 


4.86 
1.7 


140 
120 


4.99 


640 


2p:2p 


P10C18U4 


PC 


vL 


A 


25 


i 






VWvL 


A 


25 


1.2 


480 


14 


59 


17 


190 


7.6 


200 


1.7 


110 






i 






VWvL 


L 


40 


1.08 


12 


14.4 


1.4 


16.1 


8.9 


7.36 


8.7 


1.69 


4.4 






s:s 


REF. 




PPM 


L 


40 


1.08 


40 


14.6 


4.3 


16.3 


12 


7.44 


12 


1.69 


11 






lp:lp 


REF. 


drift 


VWvL 


L 


50 


1.10 

2.5 


21 
400 


20.3 
12 


1.3 

42 


52.1 
23 


68 
120 


21.2 
10 


140 
130 


3.30 
2J. 


9.7 
100 


4.99 


44 


lp:lp 


P10C18U6 


PC 


vL A 25 


i 






VWvL A 25 


2.4 


710 


18 


61 


27 


220 


12 


230 


M 


34 






i 






vL 
VWvL 


L 40 
L 40 


3.17 

2.41 


18 
460 


16.8 
20.0 


1.6 

27 


23.9 
42.6 


4.8 
190 


10.8 
19.4 


4.8 
190 


2.45 
3.02 


3.2 
120 






s:s 
2p:2p 


REF. 




PPM 


L 


40 


2.34 


550 


20.9 


34 


50.1 


230 


22.9 


230 


3.26 


150 






2p:2p 




drift 


vL 


L 


50 


2.88 


500 


24.8 


22 


43.7 


750 


23.9 


1200 


4.42 


46 


-0.286 18000 


2p:2p 


REF. 




VWvL 


L 


50 


2.28 


54 
280 


23.5 
ii 


2.8 
45 


65.5 
22 


95 
50 


27.2 
5A 


260 
49 


4.38 
5.0 


16 
38 


4.40 


150 


lp:lp 


P13C14U6 


PC 


vL 


A 


32 


i 






VWvL 


A 


31 


*A 


95 


7.5 


15 


18 


53 


4.1 


54 


53 


26 






i 






vL 


L 40 


3.86 


14 


7.40 


3.3 


16.1 


11 


3.65 


11 


4.84 


5.5 






s:s 






VWvL 


L 


40 


2.83 


6.5 


7.78 


0.76 


17.2 


3.2 


3.92 


3.2 


4.73 


0.087 






s:s 


REF. 




PPM 


L 


51 


2.66 


2.0 


7.88 


0.89 


17.3 


1.8 


3.95 


1.8 


4.72 


0.54 






s:s 




drift 


vL 


L 


55 


1.79 


190 


10.2 


16 


32.0 


650 


6.56 


860 


5.93 


180 


6.51 


190 


i:i 


REF. 




VWvL 


L 


55 


2.21 


7.1 

42 


10.2 
13 


15 

7.7 


42.8 
13 


110 

27 


8.28 


200 
34 


7.78 
Z7 


9.5 
17 


4.09 


90 


s:2p 


P13C16U4 


PC 


vL 


A 


30 


S 






VWvL 


A 


30 


M 


58 


13 


5.4 


15 


21 


5J 


22 


2.9 


12 






s 






VWvL 


L 


51 


3.57 


6.1 


14.6 


0.69 


19.1 


2.9 


6.54 


2.9 


3.11 


1.4 






s:s 


REF. 




PPM 


L 


51 


3.69 


9.9 


14.7 


1.1 


19.8 


5.1 


6.78 


5.4 


3.17 


3.4 






s:lp 




drift 


vL 


L 


56 


2.47 


240 


18.6 


10 


53.0 


310 


16.6 


560 


5.07 


45 


3.34 


320 


2p:2p 


REF. 




VWvL 


L 


56 


2.78 
4.9 


9.5 
160 


17.8 
15 


0.43 
200 


54.5 
21 


6.8 
48 


16.8 


15 
500 


4.99 
3J. 


1.4 
35 


3.47 


43 


s:lp 


P13C16U6 


PC 


vL 


A 


30 


i 






VWvL 


A 


30 


5.2 


250 


15 


23 


21 


86 


H 


86 


3J. 


61 






i 






VWvL 


L 


51 


4.84 


2.9 


15.9 


1.3 


24.5 


2.6 


8.39 


2.6 


3.56 


2.3 






s:s 


REF. 




PPM 


L 


51 


4.97 


14 


16.1 


1.4 


25.5 


2.4 


8.75 


2.4 


3.65 


2.0 






s:s 






PPM 


U 


30 


4.71 


103 


15.6 


9.3 


24.4 


28 


8.35 


27 


3.45 


74 






iP-lE 




drift 


vL 


L 


50 


4.72 


400 


20.6 


20 


55.1 


530 


22.4 


880 


5.88 


21 


2.60 


970 


2p:2p 


REF. 




VWvL 


L 


55 


4.40 


14 


19.9 


0.66 


62.5 


16 


19.4 


42 


6.11 


4.5 


3.40 


110 


s:lp 






VWvL 


U 


30 


4.01 


90 


19.7 


6.5 


56.2 


300 


19.3 


600 


5.89 


25 


3.29 


220 


lp:lp 
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Table 4. Properties temporally averaged at the outer boundary for model P13C16U6 (see Sect. 13. 3) using different numbers of gridpoints. From the left columns 
specify input properties: PC or drift type, number of gridpoints ng, advection scheme (A), and grid type (g). Remaining columns specify seven properties: the 
mass loss rate (A/), terminal velocity (noo)> degree of condensation (/cond). dust/gas density ratio {pi/p), mean grain radius (rj), drift velocity {v-q), and 
type of outflow (type). All models have been calculated with the outer boundary fixed at Text = 50 /?* . For further details see the caption of Table[3] 



drift/PC ng 
model 



10" (A/) 
[M0yr-i] 
lO^r 



(«oo> 

[kms"-"- 
10^ 



(jcond/ 

[%] 

lO'^f 



(Pd/p) 

[10-1] 
lO^r 



10" in) 

[^tm] 
lO-'^f 



[kms-1] 
lO^r 



type 



PC 700 


vL 


A 


5.3 


200 


15 


25 


22 


33 


lA 


32 


3.3 


28 






1 




VWvL 


A 


6J) 


260 


15 


22 


22 


75 


lA 


75 


3.4 


52 






i 




vL 


L 


5.49 


7.4 


15.3 


25 


22.6 


5.8 


7.72 


5.8 


3.51 


2.7 






s:s 




VWvL 


L 


5.06 


20 


16.1 


2.0 


25.4 


7.1 


8.70 


7.1 


3.66 


7.4 






s:s 




PPM 


L 


4.99 


50 


16.2 


4.3 


25.6 


8.5 


8.76 


8.7 


3.65 


15 






iP-lE 


drift 


vL 


L 


4.57 


630 


21.0 


27 


45.2 


810 


23.8 1200 


5.60 


85 


2.56 1300 


2p:2p 




VWvL 


L 


4.54 


76 
160 


19.8 
15 


4.9 
20 


61.1 
21 


130 
48 


19.4 
7A 


320 
500 


6.09 


18 
36 


3.26 


140 


lp:lp 


PC 500 


vL 


A 


4.9 


i 




VWvL 


A 


5.2 


250 


15 


23 


21 


86 


lA 


860 


M. 


61 






1 




vL 


L 


4.79 


10 


15.9 


0.46 


24.4 


2.9 


8.37 


2.9 


3.55 


2.6 






s:s 




VWvL 


L 


4.84 


2.9 


15.9 


1.3 


24.5 


2.6 


8.39 


2.6 


3.56 


2.3 






s:s 


REF. 


PPM 


L 


4.97 


14 


16.1 


1.4 


25.5 


2.4 


8.75 


2.4 


3.65 


2.0 






s:s 


drift 


vL 


L 


4.72 


400 


20.6 


20 


55.1 


530 


22.4 


880 


5.88 


21 


2.60 


970 


2p:2p 


REF. 


VWvL 


L 


4.40 


14 
63 


19.9 
15 


0.66 
2.1 


62.5 
21 


16 
18 


19.4 


42 
18 


6.11 
3.3 


4.5 
8.3 


3.40 


110 


s:lp 


PC 300 


vL 


A 


M 


i 




VWvL 


A 


6J. 


65 


15 


6.3 


21 


39 


lA 


39 


3.3 


21 






i 




vL 
VWvL 


L 
L 


6.15 

4.75 


53 
4.3 


14.6 

15.7 


5.7 
0.83 


20.1 
23.8 


3.4 
1.6 


6.84 
8.14 


34 
1.6 


3.40 
3.51 


17 
1.0 






9p:9p 

s:s 




PPM 


L 


4.76 


6.1 


15.9 


1.1 


24.4 


2.3 


8.37 


2.2 


3.56 


2.5 






lp:lP 


drift 


vL 


L 


5.06 


100 


18.9 


6.5 


58.0 


200 


18.9 


420 


6.11 


12 


4.26 


230 


2p:2p 




VWvL 


L 


4.21 


120 
30 


20.0 
13 


11 
4.0 


62.5 
16 


140 
10 


20.2 
5^ 


350 
10 


6.02 


9.5 
6.0 


3.55 


190 


2p:2p 


PC 100 


vL 


A 


5.1 


s 




VWvL 


A 


5.4 


20 


15 


1.1 


21 


2.2 


73 


2.2 


3.3 


3.5 






s 




vL 
VWvL 


L 
L 


4.96 
3.66 


13 
3.3 


14.4 
13.8 


3.8 
0.18 


15.9 
17.0 


7.1 
0.13 


5.35 
5.83 


7.3 
0.14 


2.58 
2.93 


4.9 
0.24 






s:s 
s:s 




PPM 


L 


4.16 


17 


14.7 


0.91 


20.4 


2.5 


6.97 


2.3 


3.34 


3.4 






s:s 


drift 


vL 


L 


4.15 


150 


20.1 


10 


45.2 


110 


13.5 


210 


4.88 


33 


8.66 


90 


s:s 



P10C18U6 does, with a 26 per cent larger mass loss rate calculated 
here. P13C16U4 and P13C14U6 show a mass loss rate 35 per cent 
and 28 per cent larger than earlier, while the degree of condensation 
and dust-to-gas ratios are 15 per cent smaller than before, for both 
models. As will be seen in Sect.|4]P13C14U6 is very sensitive to 
the numerical method which partly justifies its deviation. The new 
P13C16U4 develops a quasi-stationary outflow (similar to that of 
the model using a non-adaptive grid) while the structure of the pre- 
viously calculated model is more irregular. Summarising four out 
of six models show a larger mass loss rate here than before. Other- 
wise, the values are fairly similar. 



A comparison between models calculated using either the 
VWvL or PPM advection schemes and a non-adaptive grid exhibits 
a very close agreement in most cases. These two volume-weighted 
second and third order spatially accurate schemes are implemented 
in different subroutines. The good agreement of the respective re- 
sults indicates correct implementations. 



4 DISCUSSION 

In this section results are discussed which are important for a phys- 
ical interpretation of wind models. Effects of numerical origin are 
treated in Sect. |4. 11 - which can be omitted by readers not concerned 
with details of the numerical method. Changes to the physical struc- 
ture, which are obtained with the improved numerical descriptions, 
are discussed in Sect.|4.2| 



4.1 Effects of numerical modifications 

The response of a stellar wind model to the different numerical 
adjustments introduced in Sect.|2]is complex. For clarity, the influ- 
ence of each individual modification is in the following discussed 
separately. Effects of the adopted advection scheme for both drift 
and PC models are treated first in Sect. l4. l.TI Thereafter differences 
in the physical structure between models using a non-adaptive and 
adaptive grid are studied in Sect. 14.1.21 Followed by a discussion 
of the influence of the number of gridpoints in Sect. 14.1.31 
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4.1.1 Dijferences due to the accuracy of advection schemes 

In this subsection the outcome is compared in detail of models us- 
ing the three advection schemes introduced in Sect. 12.31 i.e., vL, 
VWvL, and PPM. First PC models are treated which use an adap- 
tive grid; thereafter, drift and PC models using a non-adaptive grid. 

In models that use an adaptive grid shocks are resolved lo- 
cally by densely arranged gridpoints. A more accurate advection 
scheme ought not be able to improve the outcome, unless the im- 
proved accuracy is significant also in unresolved regions. Moderate 
differences (<20per cent; increasing, as well as decreasing) are 
found between properties of models calculated using either the vL 
or VWvL advection schemes. The shape of the physical structure is 
not affected - all remain irregular Two exceptions are PI3C14U6 
and P13C16U6,ig=ioo. PI3CI4U6 is a highly numerically sensitive 
model for which all properties change significantly with the advec- 
tion scheme (by 20-32 per cent); even the terminal velocity. The 
sparsely resolved model P13C16U6,i„=i()o shows a difference of 
about 30 per cent in the dust-to-gas density ratio. It seems unneces- 
sary to make further improvements to these models as the adaptive 
grid itself is found to introduce larger eiTors (see next subsection). 

With a non-adaptive grid, which does not track shocks, shocks 
are not resolved by many gridpoints and the accuracy of the advec- 
tion scheme is always important. If a numerical scheme is robust 
models using increasingly accurate advection schemes should re- 
sult in more similar results. When comparing PC models using the 
VWvL and PPM schemes, differences in the outcome are for all 
presented properties < 6 per cent, with the exception of P10C18U6 
and P13C16U6„g=ioo- For P10C18U6 the dust-to-gas density ra- 
tio is 18 per cent larger in the PPM model compared to the VWvL 
model. For the low resolution model P13C16U6,i„=io() differences 
are 6.5 per cent (in the terminal velocity) to 20 per cent (in the de- 
gree of condensation and dust-to-gas density ratio). In these cases 
the increased accuracy is important and enhances the temperature 
sensitive dust formation. The third-order spatially accurate PPM 
models are, moreover, better at preserving the physical structure 
than VWvL models. A periodic structure - although mostly of 
small amplitude - is found far out both in the gas and dust, also with 
as few as 100 gridpoints with PPM. The same structures are more 
often smeared out when using less accurate advection schemes. 

Corresponding comparisons between PC models using the vL 
and PPM advection schemes show much larger differences in out- 
flow properties; from I.3-I10per cent. For the PI3CI6U6 models 
using different number of gridpoints (Table |4ll the same changes 
are 1.3-30 per cent; being the largest for the 7ig=100 model. The 
amount of dust is, furthermore, found to increase by 4.5-1 10 per 
cent for all models. The structure is finally found to be more sta- 
ble with increased advection accuracy - tending towards stationary 
outflows with small periodic fluctuations. These stationary outflows 
are not to be confused with stationary winds, where velocities stay 
unchanged everywhere. 

For drift models using a non-adaptive grid the increased preci- 
sion of the VWvL scheme, compared to vL, has a strong influence 
on the entire wind structure. Quantitatively, the degree of condensa- 
tion is found to increase by 1.2-50 per cent when using the VWvL 
scheme, while the drift velocity changes by 3.9-57 per cent; for 
PI0CI8U6 there is even a change of sign. The vL scheme is not 
capable of calculating the correct structure, and the structure type 
is affected for all drift models. VWvL models in every case have 
a shorter periodicity, or a stationary outflow, compared to the cor- 
responding vL model. They also, in most cases, show a smaller 
variability (seen in the f number of each property). 



It is difficult to predict which model will be more affected by 
the accuracy of the advection scheme. The final conclusion is that 
one should use a scheme as accurate as possible, but nothing less 
accurate than VWvL. In particular, for drift models there is a risk 
that the physical structure vanishes completely in numerical noise 
with a less accurate advection scheme. 



4.1.2 Comparing models using a non-adaptive/adaptive grid 

The change from an adaptive to a non-adaptive grid turns out to be 
of fundamental importance to physical properties of a stellar wind. 
Although an adaptive grid equation is able to resolve a small num- 
ber of physical features, e.g., shocks (see Sect. |2.4.2l (. it is not able 
to resolve a large number; because, there are simply not enough 
gridpoints. In addition, regions away from resolved features are 
poorly resolved. In contrast, a disadvantage with a non-adaptive 
grid is that no shock will be resolved by more than one or a few 
gridpoints (depending on the added amount of tensor viscosity and 
number of gridpoints), increasing the need for a more accurate ad- 
vection scheme. Additionally, very narrow shocks, possibly hosting 
drastically different physical conditions, can also not be studied. 

The major dichotomy between physical structures formed 
with the two different grids is: an adaptive grid mostly forms irreg- 
ular structures, while a non-adaptive grid (in thepresented cases) 
forms periodic stationary outflows, or nearly sqj. This holds for 
both drift and PC models. A comparison of the temporal out- 
flow structure of model PI0CI8U6 illustrating this point is shown 
in Fig. [T] using a non-adaptive (left panels) vs. adaptive grid 
(right panels). Furthermore, the radial structure of another model, 
PI3CI4U6, is shown in Fig.|2]using a non-adaptive grid for both the 
drift and PC model. Compared to a model using an adaptive grid, 
shocks are only strong in the inner envelope and rapidly become 
weaker moving outwards. Small periodic variations of both gas and 
dust properties, at the location of the outer boundary (originating in 
the piston), are seen for all quantities in Fig.[T](left panels). Struc- 
tural shapes of the other models using a non-adaptive grid are very 
similar to those shown in Figs.[T]&|2] 

These findings strongly indicate that an adaptive grid equation 
mostly is uncapable of tracking many weak structures, simultane- 
ously providing enough resolution in unresolved regions. At least 
not in its normally adopted setup where grid weights (for the den- 
sity and energy) are set to Wj = 1.0 (see Eq.[TJ. If these weights are 
instead replaced with Wj = 10~'^ the grid is more stiff, and struc- 
tures become periodic. The number or gridpoints across a shock in 
the inner envelope, where the density changes by about I dex, is 
with ■Wj — 1Q~'^ around 8. With Wj — IQ^^ there are about 25 grid- 
points across a density jump of about 1.3 dex, although with such 
grid weights the dust component is no longer periodic, but irregular. 
10 gridpoints are sufficient to resolve a shock of 1 dex magnitude 
(cf. Sect. 12. It . With more resolved shocks inaccuracies in unre- 
solved regions become dominating (first for the dust) and the re- 
sulting physical structure is then affected by numerical errors. This 
shows that if an adaptive grid is used to track shocks in a wind, grid 
weights should be chosen with extreme care. If grid weights are 
too large it is not possible to discriminate between physically pe- 
riodic and irregular structures (concerning models using a constant 
gas opacity, see Sect. |3.4| ). If further resolved shocks are important 



^ All models that use a non-adaptive grid and do not form stationary out- 
flows use the less accurate vL advection scheme. 
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Figure 1. Temporal evolution at the outer boundary, covering a time span of 25 P (piston periods; = 36 years). The model P10C18U6 is shown in the left 
panels using a non-adaptive giid, and in the panels on the light using an adaptive grid (with grid weights Wp^e = 1-0). Drift models are drawn with solid lines 
and the corresponding PC models with dash-dotted lines, all models use the VWvL advection scheme. From the top the panels show six properties: a) the 
terminal velocity Uac, b) mass loss rate M (logarithmic), c) degree of condensation /cond. d) dust-to-gas density ratio Pd/p (log.), e) mean grain radius rj, and 
f) drift velocity v^. Note that the models using a non-adaptive grid (left) are both periodic, while the models using an adaptive grid are irregular (right). The 
significantly larger amounts of dust and larger grains formed in the drift model are clearly seen in panels c) and e). The non-zero drift velocity in f) indicates 
well decoupled gas and dust components. Also note that the periodicity of the diift and PC models (left) differ; 1 P and 2 P, respectively. Furthermore, the 
models on the left have not settled into a periodic variation for times shorter than about 35 P. For further details, see Sect. 14.1.21 
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Figure 2. Radial structure of an instant of the model P13C14U6 using a non-adaptive grid, illustrating the full modelled domain. The drift model is shown with 
a dotted filled line and the PC model (PPM) with a solid line. From the top left panels show ten properties: a) the gas velocity u, b) gas density p (logarithmic), 
c) drift velocity vj), d) dust density pi (log.), e) nucleation rate J* (log.), f) temperature Tg, g) net growth rate r^ (log.), h) degree of condensation /cond, i) 
volume integrated growth term (cf. Sect. 14.2.21 log.), and j) average grain radius rj. All plots are drawn as a function of the stellar radius Ri, (lower axis), and 
astronomical units (au; upper axis). Each dot on the filled contours represents an individual gridpoint. Grey horizontal lines are guides. The multi-periodicity 
of this model is 2 P (drift) and 1 P (PC). For further details see Sect. l4.1.2l 



for any reason, it appears necessary to use an adaptive mesh re- 
finement instead of an adaptive grid equation. W06 presents stellar 
winds calculated using such an approach. 

In this section models calculated using a non-adaptive grid (of 
grid type L [and U]) are next compared quantitatively with the cor- 
responding models calculated using an adaptive grid. Only refer- 
ence models are considered of the non-adaptive (L) models. 

Out of six sets of model parameters three show mainly de- 
creased values on the averaged properties from PC models; using a 
non-adaptive grid. The mass loss rate of P10C16U6 is 30 per cent 
lower, while the terminal velocity, degree of condensation, dust- 
to-gas ratio, and mean grain radius are 20, 30, 30, and 10 per cent 
lower, respectively. The corresponding values for P13C14U6 are 
—40, +5, —4,-4, and — 10 per cent (with the vL advection scheme 
these values are larger: —60, —30, —20, —20, and —6 per cent). 
All quantities for P IOC 1 8U4 change by 1-10 per cent. The three re- 
maining sets of PC models show mainly increased values. The mass 
loss rate of P10C18U6 is 3 per cent lower, but the terminal velocity, 
degree of condensation, dust-to-gas ratio, and mean grain radius in- 
crease by 20, 90, 90, and 70 per cent, respectively. For P13C16U4 
the corresponding values are 10-40 per cent. P13C16U6 shows dif- 
ferent changes depending on the number of gridpoints. For ng=500 
the values are: —4, +7, +20, +20, and +20 per cent, respectively. 



while they for ng=100 are: -20, -2, 



-5, and +1 per cent. 



Consequently, not only does the structure change with the use of 
the adaptive grid equation, but average values also change. In the 
presented cases by as much as 90 per cent from a model using the 
numerical setup of previous papers in the series. 



Most models have also been calculated using a uniform non- 
adaptive grid (U). Such models do not resolve shocks well and a 
lot more artificial viscosity is necessary to widen jumps; about five 
times as much is used in the central parts of the wind, typically. 
Nevertheless, the outcome is similar to that of the other models 
using a non-adaptive logarithmic grid. Differences are smaller with 
PC models than with drift models. Results of two such (PC and 
drift) models, P13C16U6, are shown in Table [3] Their averaged 
properties deviate less from the properties of the two models using 
a logarithmic non-adaptive grid, than they do from the properties 
calculated using an adaptive grid. However, since these models are 
not capable of reproducing results to the same accuracy as models 
using a logarithmic grid they are not considered here further. 
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4.1.3 Importance of sufficient global spatial resolution 

Compared to how well shocks are resolved locally, the global res- 
olution of the full model domain also influences results notably. 
If a physical structure is well enough resolved, a stable numerical 
algorithm ought to reproduce results with an increased number of 
gridpoints. Exactly what happens with fewer gridpoints depends on 
how important narrow structures - that are then unresolved - are to 
the solution. In this section the results presented in Table|4]are dis- 
cussed. To simplify the interpretation, presented PC and drift mod- 
els are compared with the respective P13C16U6n„=50() reference 
model. Models using a non-adaptive grid are discussed first. 

Shocks are always resolved with at least one gridpoint on a 
shock in the models with rig = 500 and 700; see Fig. [3] (dot- 
ted filled line) for an illustration. All properties of the drift model 
with rig = 700 are found to differ by < 4.1 per cent, and of the 
PC model by < 0.40 per cent, from the corresponding reference 
model. The larger differences between the drift models are likely 
due to the lower accuracy of the VWvL advection scheme. Changes 
to the physical structure are negligible, albeit "fluctuations" are 
slightly larger for the model using ng=700 (hence the difference 
in the characterised type in Table |4j. With rig = 300 shocks are 
not well resolved, but average values are still satisfactorily repro- 
duced, both for the PC model - < 4.3 per cent (VWvL: < 7.0 per 
cent) - and the drift model; < 4.4 per cent. Furthermore, while the 
PC model shows small periodic variations, the periodicity of the 
drift model changes from IP to 2 P. With jig = 100 differences 
are more drastic and the PC model shows decreased values in all 
quantities (< 20 per cent; VWvL: < 33 per cent). How unresolved 
shocks with Ug = 100 appear is illustrated with the drift model 
P13CI6U6ng=ioo in Figs.[3ji & c; compare the solid line with the 
dotted line, which uses rig = 100. The drift velocity reaches higher 
values in this model (40 km s^^, and above), than it does when us- 
ing more gridpoints and a more accurate advection scheme (than 
vL), see Fig. |3}) (2 i?* ^ r < 4 R^^^, compare with the dotted line 
using ng=700). 

Comparing the PC models using the original vL advection 
scheme and an adaptive grid with the reference PC model these 
all show an increased mass loss rate and decreased numbers on the 
remaining quantities. The differences are the smallest for rig = 700 
(5.0-14per cent), intermediate for rig = 300 (10-21 per cent), and 
the largest for Ug = 100 (2.6-37 per cent). The latter model com- 
pared to the other adaptive grid models also forms a stationary out- 
flow - since there are too few gridpoints available to resolve shocks. 



4.2 Physical properties of a dust-driven stellar wind 

Previous sections have focused on the implementation and analysis 
of numerical improvements. In this section physical structures of 
stellar winds are studied, which form using these improvements. 
An emphasis is put on effects of gas-dust drift. 

Compared to (non-drift) models that assume position coupling 
(PC), drift models have been found to show many differences. Ef- 
fects of d rift, as found previously, can be summarised as follows 
(sect. 4 in lSandinluOOSi , contains a more detailed summary): 

• Larger variations occur in the physical structure, in particular 
in the dust component. Dust accumulates to locations of shocks, 
where densities are larger; and the gas-dust interaction stronger. 

• Dust and gas components are due to strong shocks, also 
present in the outer envelope, well coupled. The drift velocity typ- 
ically assumes values in the range 0.1 ^ wd ^ 40 kms~^; larger 
values occur between shocks where the density is lower. 




Figure 3. Radial structure of an instant of the inner wind formation re- 
gion of the two drift models P13C16U6nj=700 (dotted filled line) and 
P13C16U6n„=ioo (diamonds, o); dots and diamonds show positions of in- 
dividual gridpoints. Three quantities are shown: a) the gas velocity u, b) 
drift velocity vj), and c) gas density p (log.). Note the wide (unresolved) 
shocks in the model using rig = 100. For further details see Sect. 14. 1. 31 



• Larger amounts of dust form when assuming drift-dependent 
dust formation; specifically grain growth is enhanced. 

• Mass loss rates are in many cases found to be lower. It is as- 
suming drift harder to form certain wind models, which are found 
to form weak winds when assuming PC; weak winds are defined to 
be those where the terminal velocity Uoa < 10 kms^^. 

Moreover, most of the Planck-mean and constant-opacity (drift and 
PC) models are found to form an irregular (instead of a periodic) 
structure; although constant-opacity models are occasionally peri- 
odic. Due to larger variations in the dust component drift model 
structures are more irregular than those of PC models. Conclusions 
are sometimes ambiguous. For example, with an irregular structure 
it is difficult to see to what extent dust accumulates to shocks, and 
occasionally a drift model is found to generate a higher mass loss 
rate. In all cases, however, drift models form more dust. 

Disentangling the interplay between a large number of phys- 
ical processes (and a numerical dependence) is difficult with a 
strongly irregular structure. Differences are with the new more peri- 
odic models presented here found to be more clear and pronounced. 

Wind structures, variation patterns, and amount of dust formed 
are, furthermore, highly dependent on the adopted physical as- 
sumptions and parameters describing the wind model. Absolute 
numbers are expected to change with the input. Time-dependency, 
opacities, piston properties, dust properties, and stellar parameters 
are such important assumptions. 

In the following discussion only models using the best avail- 
able advection scheme and a non-adaptive grid are used; these mod- 
els are marked with "REF." in Table [3] Issues of wind variabil- 
ity, drift velocity, and gas-dust decoupling are discussed next in 
Sect. 14.2.11 Thereafter, revised findings on the efficiency of dust 
formation are treated in Sect. |4.2.2l followed by a quantitative com- 
parison of how remaining averaged properties differ between drift 
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and PC models in Sect. 14.2.31 Finally, some issues for model im- 
provement are discussed in Sect. 14.2.41 

4.2. 1 On the drift velocity and variability of model structures 



The critical finding of Sect. l4.1.2l was that all (revised) wind struc- 
tures are periodic. On top of such periodic structures the gas com- 
ponent in every case shows relatively small variations. In order to 
see this, compare the gas velocity and density structures in Figs. [2^ 
& b and Figs.|4^ & d with, e.g., figs. 2a & b in Paper III (which illus- 
trates a model with an irregular structure). Except in the wind accel- 
eration region and below (i.e., for r > 3 i?*), gas shocks are weak - 
typically the change across a shock is < 0.7 dex in the density; this 
value is even smaller in PC models - and rapidly decrease in am- 
plitude moving outwards. The dust component likewise shows very 
small variations in PC models; an exception is P10C18U6, which 
shows large variations for both components (and also strong dust 
formation). 

Drift models show a much more variable dust component than 
PC models do, confirming the previous finding of variability. Un- 
like earlier, the dust component is now much more decoupled from 
the gas. A consequence partly linked with shocks, which are weaker 
than earlier when structures were irregular; the gas-dust interaction 
is therefore also weaker. Although dust (primarily) forms in the 
dense environment of a shock, it then drifts ahead of the gas (since 
vd > 0). This partial decoupling is shown for model P13C16U6 
in Fig. |4] Comparing the locations of the two dust "shells" in the 
degree of condensation (Fig.|4};) with the locations of peaks in the 
gas density (Fig.|4ji; and velocity, panel a) for 16 -R* ^ r ^ 25 i?*, 
it is seen that they do not match. Also note that the drift veloc- 
ity (Fig.|4j)) in the same region attains peak values in front of the 
(weak) gas shocks - sometimes displaced from the locations of dust 
shells. Differing locations of maxima in various gas and dust prop- 
erties indicate that an interaction between the two components is 
still taking place. An interaction "breaking" the dust relative to the 
gas (this is seen in the decreased drift velocity with radius. Fig. [2};). 

What is not seen in the presented figures is that the entire drift 
velocity structure is oscillating at a large amplitude, on top of the 
variations with radius seen in, e.g., Fig.|4j). At the outer boundary 
the oscillations look like in Fig. [if (left panel); also see the super- 
posed periodic variations in Paper II (fig. 5 [bottom right panel]). 
Such oscillations demonstrate how tightly dust is bound to the ra- 
diation field. 

Another new result is that the dust component in the outer en- 
velope loses its "patterned" structure (i.e., dust shells are smeared 
out spatially); as is seen in the degree of condensation in Fig. [2jr 
(for r > 30 7?*). The same behaviour is found in the other drift 
models as well. With current numerical limitations it is difficult 
to model the region beyond about r = 50 i?* due to the locally 
rapidly decreasing grid resolution of the logarithmic grid. Hence, 
any statement on what happens to the variability pattern further out 
is uncertain. The same dust pattern is also diffused away when us- 
ing a uniform grid, but there the grid resolution problem is even 
more alarming; but then in the inner envelope. It consequently ap- 
pears that also the dust component of drift models ultimately form 
stationary outflows, at larger distances from the star. 

The drift velocity is for a vast majority of circumstances (mea- 
sured spatially and temporally) found to achieve moderate val- 
ues; wd < 10 kms^^. Temporally averaged values of the drift 
velocity are shown as a function of mass loss rate in Fig. |5}l; re- 
vealing a trend of slightly decreasing values with increasing mass 
loss rate. The highest values occur in the wind acceleration re- 



gion 3 -R* < r < 10 -R*, decreasing outwards. Non-thermal sputter- 
ing, requiring drift velocities of about 40 km s~^, does not play a 
role (a statement likely to change if a grain size distribution is used; 
cf. Sect. l4.2.4t . Moreover, the primary molecule of importance to 
the grain growth process is, under these circumstances, C2H2 (see, 
e.g.. Paper III). Never is the drift velocity subsonic outwards of the 
dust forming region. The speed ratio (Sd\ e.g.. Paper I, eq. 11), is 
found to take on values of So ~ 2-3 (4 for P10C16U6) in most 
parts of the envelope. Referring to Paper I (fig. 1) it is found that 
the relative error of the used drag coefficient {Cif^; ibid., eq. 23) 
is about 0.5-1.5 per cent for these values. With the high-velocity 
approximation (C^^; ibid., eq. 24) corresponding errors would be 
10-20 per cent, which is why its use is not recommended. 

4.2.2 On the amount of formed dust 

Dust formation can occur where the density is high enough to al- 
low for abundant collisions between gas (and dust) particles. A 
simultaneous requirement is that the temperature is low enough 
that dust grains do not evaporate. The resulting region where grain 
growth is efficient has been found to be relativ ely small - about 
2 Rj, 5, '' 5, 5 Ri, (in the inner en velope; e.g., lOail & Sedlmavn 
ll987l ; lDomini"3ll992l ; lH6fneill2007h . 

Using a stationary wind model, that included drift-dependent 
dust formation (and a grain size distribution), KS97 (see sect. 4.1) 
found, in one of two presented models, that the dust-to-gas mass 
loss ratio ({Mi) / {M}), increased by about 50 per cent due to drift. 
Such a "flux" ratio is more appropriate to use than the degree of 
condensation when comparing amounts of formed dust, since the 
dust component is dynamically diluted compared t o the gas with a 
non-zero drift velocity (also see lKrtiger et aljl994h . 

For an averaged grain size distribution the dust-to-gas mass 
loss ratio is. 

Mi _ pdU + VD 

M p u ' 

where it is the gas velocity. Resulting ratios for all presented mod- 
els are shown in Fig.[5^ as a function of the mass loss rate (M). As 
the same figure (lower panel) illustrates the increased amounts of 
dust due to drift, 160-250 per cent, which are fed to the interstellar 
medium, are larger than the corresponding increase in the degree 
of condensation (1 10-220 per cent); see Fig.jSj) (lower panel). Dif- 
ferences between mass loss ratios and dust-to-gas density ratios, 
due to dilution, are about 17-40 per cent. One model, P10C18U6, 
shows an efficient dust formation in the PC model (partly due to a 
high piston velocity and carbon-to-oxygen ratio), which is why the 
increase to the drift model (41 per cent) is smaller compared to the 
other five models. Grain radii are, simultaneously, 34-95 per cent 
larger (Fig. |5j:). 

Comparing the mass loss ratios of the bulk of the models 
(2.6-3.5) with the value of KS97 ((Md)/{M)«1.5), this differ- 
ence is to be expected since the adopted stellar parameters are very 
different. Dynamical effects are less well accounted for in non- 
variable velocity structures of stationary winds. Features such as 
shocks, which promote dust formation, cannot form. 

The negative slope of the dust-to-gas mass loss ratio with in- 
creasing mass loss rate seen in Fig. [5^ (lower panel) seems to indi- 
cate that dust formation (and therefore the envelope) is less affected 
by drift at larger mass loss rates (and hence larger densities); this 
is also noticed by KS97. Under such circumstances the gas-dust 
coupling is stronger, which leads to smaller drift velocities and sup- 
posedly more similar drift and PC models. On the contrary, at lower 
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Figure 4. Radial structure of the model P13C16U6 using a non-adaptive grid, illustrating the inner model domain. The drift model is shown with a dotted 
filled line and the PC model (PPM) with a solid line. From the top the panels show four properties: a) the gas velocity u, b) drift velocity I'd, c) degree of 
condensation fcoad, and d) gas density p (logarithmic). All plots are drawn as a function of the stellai' radius R^, (lower axis), with a complementing measure 
in astronomical units (au; upper axis). Dots on the filled contours represent individual giidpoints. Note that the gas and dust components are loosely coupled, 
compare the locations of maxima in the gas density, drift velocity and degree of condensation for 17 i?* sCr ^25/?*. For further details see Sect. l4.2."i1 



densities the coupling is weaker. Too few models have been calcu- 
lated to permit a reliable estimate of the model behaviour outside 
the current sets of model parameters. It is likely that the range of 
mass loss rates where drift (and dust) is important is limited. This 
statement could be different if it turns out that models with higher 
mass loss rates form irregular winds (e.g., by the use of other opac- 
ities, cf Sect. 14.2.4b . 



Time scales for grain growth are longer in the outer envelope 
than in the inner, since densities there are lower. The volume of 
a shell in a radial interval increases with radius in a spherical ge- 
ometry, which is why the volume-integrated amounts of formed 
grains still could be significant. In Fig. [2j the volume integrated 
dust growth (source) term of the moment equation corresponding 
to the dust (number) density is shown (i.e., the first source term 
on the right hand side of the K3 moment equation; see, e.g.. Pa- 
per III, eq. 2). This term is a lot flatter than the net growth rate 
(tq ^ , Fig.|2^), and decreases by less than a factor of two across the 
modelled envelope. Comparing the peaks of formation in the drift 
model the same decrease is even smaller; about a factor of one. In 
Fig. |2^ the growth rate of the drift model is about 0.7 dex larger 
than it is in the PC model for r > 10 i?*. The grain growth source 
term can be integrated across the envelope for both the drift and PC 
model in the shown instance. Doing this it is found that 12 per cent 
(PC: 8.6 per cent) of the total amount of formed dust forms at radii 
r ^ 10 7?*. It should be noted, however, that the dust formation 
efficiency changes during a pulsation period; although these ratios 
stay about the same for this model. 



4.2.3 On the change of averaged properties 

The mass loss rate is in every presented case, but one, found to be 
slightly lower in drift models; the reason being a less than complete 
gas-dust coupling. Compared to the PC models the change is about 
— 1.9-27 per cent. Another new result is that the terminal velocity 
for all drift models is larger than in the corresponding PC model; by 
about 12^5 per cent in the presented cases (Fig. [5^). The combina- 
tion of a lower mass loss rate and larger terminal velocity indicates 
that the gas in a drift model is more tenuous than in a PC model 
(compare the drift and PC model gas densities in Fig.[2j3). 

The range of increased values of the terminal velocity seems 
to be in accord ance with observed velocity distributions (see, e.g., 
IOlofssonll2003l . sect. 7.6.2 and references therein). These tend to 
peak towards higher values for (dust-enshrouded) IR C-Stars. 



4.2.4 Points for improvement with more detailed physics 

Drift models are now accurate enough that conclusions of topics 
first approached in Paper I would benefit of a new study. Such top- 
ics include the role of specular vs. diffusive collisions (in the mo- 
mentum and energy transfer) and the influence of different expres- 
sions of the drag force. Current models also leave a few additional 
issues open for further (and more difficult) improvements. Three 
such improvements are immediately apparent: using a grain size 
distribution in place of a mean, switching to frequency-dependent 
opacities, and using a two- or three-dimensional approach in place 
of the current one-dimensional. How these improvements could af- 
fect results are briefly discussed next. 

Two studies, using stationary model formulations, have been 
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Figure S. Illustration showing how properties (from Table [5) of drift and PC (reference) models compare. Upper panels show absolute numbers, and lower 
panels the drift-to-PC ratio of the respective property (given in per cent). From the left five properties are shown: a) terminal velocity (moo), b) degree 
of condensation {fcond), c) mean grain radius (rj), d) drift velocity (vj^), and e) dust-to-gas density ratio (pd/p)(left axis) & dust-to-gas mass loss ratio 
(Ali) / (M) (right axis); all drawn as a function of the mass loss rate (M). Error bars in the upper panels show weighted relative fluctuation amplitudes (r). 
Moreover, drift (PC) model values are shown with dark grey (grey) bullets, • ( ). Values of drift and PC models using the same stellar parameters are connected 
with grey lines. In panel e) black bullets (•) denote mass loss rate ratios and © indicates the ratio found by KS97 for their low mass loss model. Note that all 
presented models, with one exception, form drastically larger amounts of dust when assuming drift. For further details see Sect. l4.2.2l 



carried out earli er where the influenc e of a grain size distribution 
has been treated JDominik et al.l 19891 : KS97). Since so much larger 
amounts of dust can form with current time dependent drift models, 
it is of high interest to see what the outcome is when appropriately 
considering a grain size distribution; instead of a mean. In reality 
grains of different size are affected by the radiation pressure from 
the central star differently. Large grains have a large cross section, 
absorb more radiation, and move at higher velocities. They could 
drift so fast that further growth is inhibited by non-thermal sputter- 
ing (see, e.g.. Paper III, sect. 2.2.2). Smaller grains drift slower and 
thereby have more time to grow efficiently. Forming a grain size 
distribution similar to that found in pre-solar meteorites ( see, e.g., 
iBematowicz et alj2005l;lNuth et alj2006l ; lAndersei]|2007l and ref- 
erences therein), it appears that drift could be an important part of 
the explanation. 



As was emphasised in Sects. 1X21 & l3.4l opacities strongly in- 
fluence properties of the wind structure. It is important to study 
how results of this paper - that use grey molecular opacities - com- 
pare to result s calculated using f requency-dependent opacities (as 
introduced bv lHofner et al.ll2003h : it is difficult to predict what the 
outcome would be. If it is found that such models are also periodic, 
permitting well decoupled gas and dust components, and showing 
fluctuations of moderate amplitude, qualitative results should show 
some agreement. 

Larger variations and patterned structures are present in the 
dust component, not the gas. From the spatial distribution of dust 
shells in the ID drift models presented here - where gas is smoothly 
distributed - it seems highly likely that the dust component would 
also be more affected than the gas in 2D and 3D models. Studies 



like those performed bv lFrevtag & Hofneii ( 120030 and W06, where 
the formation of clumps can be studied, are required to show how. 
Those models do not, however, include any formulation for gas- 
dust drift currently, and will unlikely be able to reproduce an as 
varying dust component as found here. In this context it is worth 
emphasising that clumps and structures studied by W06 (in 2D) 
are calculated using a constant gas opacity, for which outflows in a 
majority of cases are found to be irregular (see Sect. 13.41 and ibid, 
sect 3.1). If the understanding of how model structures form is to 
improve significantly it seems necessary to use 2D or 3D models, 
realistic opacities, and drift simultaneously. 



5 CONCLUSIONS 

In recent years models of dust-driven winds from AGB stars have 
been made increasingly realistic with the refinement and addition of 
more physics. These models must handle many complex non-linear 
and time-dependent processes that describe a combined evolution 
of gas, dust, and a strong radiation field. Solving the resulting ra- 
diation hydrodynamic system of equations puts large demands on 
the numerical method used, which has to correctly treat, e.g., nar- 
row shocks, and time scales covering many orders of magnitude. 
Results of previous studies (using the same modelling approach as 
in earlier papers of this series) show a majority of wind structures 
that are irregular (i.e., mildly chaotic). There is no thorough study 
where the occurrence and extent of such irregularities are clarified. 
The purpose of this paper was to fill this gap by studying the 
influence of the numerical meth od on wind models that use the 
so-called adaptive grid equation l IDorfi & Drurvll 19870 . The study 
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was bas ed on th e model description introduced in lSandin & Hofnen 
( l2003iJlbl |2004 Paper I-Paper III). Clarifying the numerical in- 
fluence all physical assumptions and parameters were kept un- 
changed. As before the emphasis was on improving the understand- 
ing of effects of gas-dust drift. In order to reassess previous results 
models were calculated both allowing drift (drift models), and not 
allowing drift (position coupled [PC] models). 

An important result of this study is that an adaptive grid equa- 
tion is not in general capable of tracking a large number of shocks 
without introducing numerical errors into the solution. All pre- 
sented and revised models using a non-adaptive grid (which does 
not track shocks) were found to form periodic structures with a 
stationary outflow showing small amplitude fluctuations (originat- 
ing in the piston). This was found to be the case both for drift and 
PC models. The same structures were all irregular when the grid 
was instead adaptive (tracking shocks). Although shocks are well 
resolved with an adaptive grid, regions between shocks simultane- 
ously become unadequately resolved; because, there is only a fixed 
amount of gridpoints. In drift models dust is not bound to the gas 
and might require an appropriate resolution also between shocks. 
Using such a powerful tool as the adaptive grid equation without 
first assuring that its numerical influence does not dominate over 
the physical structure it is impossible to tell whether a modelled 
wind is physically irregular, or if variations are of numerical origin. 

In attaining a high accuracy it is, moreover, necessary to 
use a precise advection scheme and a sufficient number of grid- 
points. Detailed studies showed that smaller amounts of dust 
form with unresolved shocks. Incorrect structures also formed 
with a less accurate advect ion scheme. Using the PPM scheme 
JColella & Woodward 1 1984) . numerical errors were found to be 
small; at the per cent level. Compared to previously calculated 
models, differences in individual model properties were overall as 
large as 100 per cent. 

It is important to point out that absolute values of results in this 
study are highly dependent on adopted physical assumptions and 
parameters, as well as the range of stellar parameters used. Qual- 
itatively, the results are comparable with, or applicable to, mod- 
els that also turn out to form periodic structures. An application to 
models which instead form irregular structures is not meaningful, 
but requires a separate study. 

Many physical results of previous papers in this series were 
confirmed, but there were also new findings. The gas compo- 
nent typically showed small variations. Unlike earlier there were 
no shocks in the outer envelope. The dust component, however, 
showed large variations in drift models, confirming previous re- 
sults. Although the periodic variation pattern diffused away at 
larger radii; and it consequently appears that also the outflow of 
the dust component becomes stationary. Moreover, unlike previ- 
ous results the dust component was less tightly coupled to the gas. 
Typically the drift velocity was fairly low in the entire envelope, 
^D J; 10 kms^^. The highest values were assumed in the inner 
envelope, decreasing outwards to 3 < «d ;$ 5 km s~^ at the outer 
boundary (at about 50 stellar radii). 

Drastically larger amounts of dust formed in drift models com- 
pared to PC models; the increase was 160-250 per cent in every 
case but one (for which the increase was 40 per cent). It was, fur- 
thermore, found that mass loss rates typically were lower in drift 
models, by < 30 per cent. Terminal velocities were, to the contrary, 
larger, by 10-50 per cent. Winds of drift models were consequently 
more tenuous than those of PC models. Altogether results (again) 
show that neglecting gas-dust drift creates imprecise outcome. 

The modelled winds were neither the least, nor the most mas- 



sive. Partly this is a consequence of the Planck mean absorption 
coefficients used, which result in less dense wind structures, than 
would be the case had instead frequency-dependent opacities been 
used. Partly because weak winds are difficult to calculate with drift. 
It is highly likely that the range of mass loss rates where drift (and 
dust) is important is limited. Finally, there are issues for model im- 
provement where a closer study would make the understanding of 
the wind formation mechanism more complete. Three such issues 
addressed were: using a grain size distribution instead of a mean 
size, using frequency-dependent opacities, and the role of a two- or 
three-dimensional modelling approach. 
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